<!DOCTYPE html>
<html lang="zh-cn">
<head>
    <meta charset="utf-8" />
    <title>椭圆函数, 超几何函数</title>
    <link rel="stylesheet" type="text/css" href="../../css/note.css" />
</head>
<body>

<h2>椭圆函数</h2>

<p class="remark">
    作换元 `t = pi/2 - x`, 显然有
    <span class="formula">
        `int_0^(pi/2) f(sin t) dt = int_0^(pi/2) f(cos x) dx`,
    </span>
    这一结论在下文直接应用, 不再说明.
</p>

<h3>第一类椭圆积分</h3>

<div class="img sm" title="K(k)">
  <canvas id="plt-elliptic" width="150" height="150"></canvas>
</div>

[<a href="https://zhuanlan.zhihu.com/p/130428825">原文来自知乎</a>]

<p class="definition">
	<b>第一类椭圆积分</b>
  完全椭圆积分:
	<span class="formula">
		`K(k) := int_0^(pi/2) ("d"theta)/sqrt(1-k^2 sin^2 theta)`,
		`quad "Re"(k) lt 1`.
    <span class="label">(第一形式)</span>
	</span>
  其中变元 `k` 称为<b>模</b>.
  令 `t = sin theta` 得到:
  <span class="formula">
    `K(k) = int_0^1 dt/(sqrt(1-t^2) sqrt(1-k^2 t^2))`.
    <span class="label">(Jacobi 形式)</span>
  </span>
  利用下文的 Landen 变换得到:
	<span class="formula">
		`K(z) = int_0^(pi/2) ("d"theta)/sqrt(1+z^2 +- 2z cos 2 theta)`.
    <span class="label">(第二形式)</span>
	</span>
  不完全椭圆积分:
  <span class="formula">
    `F(varphi, k) := int_0^varphi ("d"theta)/sqrt(1-k^2 sin^2 theta)`.
  </span>
</p>

<p class="proposition">
	<b>Euler 变换</b>
	<span class="formula">
		`K("i"x) = 1/sqrt(1+x^2) K(sqrt(x^2/(1+x^2)))`.
	</span>
</p>

<p class="proof">
	<span class="formula">
    右边
		`= 1/sqrt(1+x^2) int_0^(pi/2) ("d"theta)/sqrt(1-x^2/(1+x^2) cos^2
		theta)`
		`= int_0^(pi/2) ("d"theta)/sqrt(1+x^2-x^2 cos^2 theta)`
		`= int_0^(pi/2) ("d"theta)/sqrt(1+x^2 sin^2 theta)`
    `=` 左边.
	</span>
</p>

<p class="proposition">
  <b>Landen 变换</b>
  取共轭模 `k' = sqrt(1-k^2)`.
  令 `z = (1-k')/(1+k')`, 则
  <span class="formula">
    `K(z) = (1+k')/2 K(k)`
    `= 1/(1+z) K((2sqrt z)/(1+z))`.
  </span>
</p>

<p class="proof">
  令
	<span class="formula">
		`t = (1+k') (s sqrt(1-s^2))/sqrt(1-k^2 s^2)`,
	</span>
  求微分,
	<span class="formula">
		`dt = (1+k')/sqrt(1-s^2)
		(1-2s^2 + k^2 s^4)/(1-k^2 s^2)^(3//2) "d"s`,
	</span>
  从 Jacobi 形式出发, 令 `z = (1-k')/(1+k')`,<br>
  <span class="collapse">
    由
    <span class="formula">
      `1-k^2 s^2 - (1+-k')^2 s^2 (1-s^2)`
      `= 1 - k^2 s^2 + (s^4 - s^2) (1 + {:k':}^2 +- 2 k')`
      `= 1 - 2(1+-k') s^2 + (1+-k')^2 s^4`
      `= [1 - (1+-k') s^2]^2`,
    </span>
    有
    <span class="formula">
      `sqrt(1-t^2) = (1 - (1+k') s^2)/sqrt(1-k^2 s^2)`,
      `quad sqrt(1-z^2 t^2) = (1 - (1-k') s^2)/sqrt(1-k^2 s^2)`.
    </span>
    当 `t = 0` 时, 相应的 `s = 0` 或 `1`;
    当 `t = 1` 时, 解二次方程知, 相应的 `s = 1/sqrt(1+k')`.
  </span>
	<span class="formula">
    `K(z)`
		`= int_0^1 dt/(sqrt(1-t^2) sqrt(1-z^2 t^2))`
		`= (1+k') int_0^(1//sqrt(1+k')) ("d"s)/(sqrt(1-s^2) sqrt(1-k^2 s^2))`
		`= (1+k') int_(1//sqrt(1+k'))^1 ("d"s)/(sqrt(1-s^2) sqrt(1-k^2 s^2))`
		`= (1+k')/2 int_0^1 ("d"s)/(sqrt(1-s^2) sqrt(1-k^2 s^2))`.
	</span>
</p>

<p class="proposition">
    定义
	<span class="formula">
		`M(a, b) := int_0^(pi/2)("d"theta)/sqrt(a^2 cos^2 theta+b^2 sin^2 theta)`,
	</span>
	则
	<span class="formula">
		`M((a+b)/2, sqrt(a b)) = 2/(a+b) K((a-b)/(a+b)) = M(a, b)`.
	</span>
    取 `a, b = 1 +- k`:
    <span class="formula">
        `M(1, sqrt(1-k^2)) = K(k) = M(1+k, 1-k)`.
    </span>
</p>

<p class="proof">
	首先
	<span class="formula">
		`a^2 cos^2 theta + b^2 sin^2 theta`
		`= a^2 - (a^2-b^2) sin^2 theta`,
	</span>
	故
	<span class="formula">
		`((a+b)/2)^2 cos^2 theta + a b sin^2 theta`
		`= ((a+b)/2)^2 - ((a-b)/2)^2 sin^2 theta`.
	</span>
    于是
	<span class="formula">
		`M((a+b)/2, sqrt(a b))`
		`= int_0^(pi/2) ("d"theta) / sqrt(((a+b)/2)^2 - ((a-b)/2)^2 sin^2 theta)`
        `= 2/(a+b) K((a-b)/(a+b))`.
	</span>
	这就证明了第一个等号. 接下来利用 `K(z)` 的第二形式,
	<span class="formula">
		`2/(a+b) K((a-b)/(a+b))`
		`= 1/(a+b) int_0^pi [1 + ((a-b)/(a+b))^2 - 2(a-b)/(a+b) cos
		theta]^(-1/2) "d"theta`
		`= int_0^pi [(a+b)^2 + (a-b)^2 - 2(a-b)(a+b) cos theta]^(-1/2)
		"d"theta`
		`= int_0^pi [a^2(1-cos theta) + b^2(1+cos theta)]^(-1/2) ("d"theta)/sqrt 2`
		`= int_0^pi (a^2 sin^2{:theta/2:} + b^2 cos^2{:theta/2:})^(-1/2)
		("d"theta)/2`
		`= M(a, b)`.
	</span>
	这就证明了第二个等号.
</p>

<p class="theorem">
	<b>(Gauss, 1799)</b>
	设数列 `{a_n}`, `{b_n}` 满足: `0 lt b_0 lt a_0`, 且
	<span class="formula">
		`a_(n+1) = (a_n + b_n)/2`, `quad b_(n+1) = sqrt(a_n b_n)`,
		`quad n ge 0`,
	</span>
	则两数列收敛到同一极限. 此极限称为
	`a_0` 与 `b_0` 的<b>算术几何平均值</b> `"agm"(a_0, b_0)`, 我们有
	<span class="formula">
		`"agm"(a_0, b_0) M(a_0, b_0) = pi/2`.
	</span>
    从而得到高效求 `K(k)` 的公式
    <span class="formula">
        `K(k) = pi/2 1/("agm"(1, sqrt(1-k^2)))`.
    </span>
</p>

<ol class="proof">
	<li>先证极限存在.
		由均值不等式,
		<span class="formula">
			`b_n lt sqrt(a_n b_n) lt (a_n+b_n)/2 lt a_n`,
		</span>
		故 `{[b_n, a_n]}` 形成一组闭区间套, 且
		<span class="formula">
			`|a_(n+1) - b_(n+1)| lt 1/2|a_n - b_n|`.
		</span>
		因此闭区间套的长度趋于零. 由闭区间套定理知, 两数列的极限为
		<span class="formula">
			`A = nnn_n [b_n, a_n]`.
		</span>
	</li>
	<li>注意椭圆积分 `M(a, b)` 的值在此数列的迭代下保持不变, 有
		<span class="formula">
			`M(a_0, b_0) = M(a_1, b_1) = cdots = M(A, A) = pi/(2 A)`.
		</span>
	</li>
</ol>

<p class="proposition">
	<b>Gauss 常数</b>
	<span class="formula">
		`G = "agm"(1, sqrt 2)^-1 = (Gamma^2(1/4))/(2pi)^(3/2)`.
	</span>
</p>

<p class="proof">
  `pi/2 G`
  `= pi/2 "agm"(1, sqrt 2)^-1`
  `= K("i")`
  `= int_0^1 dt/(sqrt(1-t^2) sqrt(1+t^2))`
  `= 1/4 B(1/4, 1/2)`
  `= (sqrt pi Gamma(1/4)^2 sin{:pi/4:})/(4 pi)`.
</p>

<p>一道经典例题: <a href="../../physics/mechanics/1.html">`K(z)` 与单摆周期</a>.</p>

<h3>第二类椭圆积分</h3>

<p class="definition">
    <b>第二类椭圆积分</b>
    <span class="formula">
        `E(z) = int_0^(pi/2) sqrt(1-z^2 sin^2 theta) "d"theta`
        `= int_0^1 sqrt((1-z^2 t^2)/(1-t^2)) dt`.
    </span>
</p>

<p class="corollary">
    <b>Euler 变换</b>
    `E("i"x) = sqrt(1+x^2) E(sqrt(x^2/(1+x^2)))`.
</p>

<p class="example">
    <b>椭圆的周长</b>
    记椭圆的离心率为 `k = sqrt(a^2-b^2) // a`, 则椭圆周长
    <span class="formula">
        `L = 4 int_0^(pi/2) sqrt(a^2 cos^2 theta + b^2 sin^2 theta)
        "d"theta`
        `= 4 a int_0^(pi/2) sqrt (1 - k^2 sin^2 theta) "d"theta`
        `= 4 a E(k)`.
    </span>
</p>

<p class="example">
    求 `int_0^1 sqrt((1+x^2)/(1-x^2)) dx`.
</p>

<p class="solution">
    原式等于
    <span class="formula">
        `int_0^1 (1+x^2)/sqrt(1-x^4) dx`
        `= int_0^1 1/sqrt(1-x^4) dx + int_0^1 x^2/sqrt(1-x^4) dx`
        `= I_1 + I_2`.
    </span>
    令 `t = x^4`, `dt = 4x^3 dx`, `dx = dt/(4 t^(3//4))`,
    <span class="formula">
        `I_1 = 1/4 int_0^1 t^(-3//4)(1-t)^(-1//2) dt = 1/4 B(1/4, 1/2)`.
    </span>
    类似 `I_2 = 1/4 B(3/4, 1/2)`.<br/>
    另一方面, 令 `x = sin theta`, 原式等于
    <span class="formula">
        `int_0^(pi/2) sqrt(1+sin^2 theta) "d"theta`
        `= E("i") = sqrt 2 E(1/sqrt 2)`.
    </span>
</p>

<p class="proposition">
  <b>`E` 和 `K` 的导数</b>
  [来自 菲赫金哥尔茨 第二册511目]
  <span class="formula">
    `("d"E)/("d"k) = (E-K)/k`,
    `quad ("d"K)/("d"k) = (E//{:k':}^2 - K)/k`.
  </span>
  其中共轭模 `k'` 满足 `k^2 + {:k':}^2 = 1`.
  进一步可得 `E` 满足的微分方程
  <span class="formula">
    `("d"^2 E)/{:"d"k:}^2 + 1/k ("d"E)/("d"k) + E/{:k':}^2 = 0`.
  </span>
</p>

<p class="proof">
  直接计算:
  <span class="formula">
    `("d"E)/("d"k) = "d"/("d"k) int_0^(pi/2) sqrt(1-k^2 sin^2 t) dt`
    `= int_0^(pi/2) (-k sin^2 t)/sqrt(1-k^2 sin^2 t) dt`
    `= (E-K)/k`.
  </span>
  `K` 的导数则不太简单: 先验证等式
  <span class="formula">
    `"d"/dt (sin t cos t)/sqrt(1-k^2 sin^2 t)`
    `= 1/k^2 [(1-k^2 sin^2 t)^(1/2) - (1-k^2)/k^2 (1-k^2 sin^2 t)^(-3/2)]`,
  </span>
  于是
  <span class="formula">
    `int_0^(pi/2) (1-k^2 sin^2 t)^(-3/2) dt`
    `= E//{:k':}^2`.
  </span>
  从而
  <span class="formula">
    `("d"K)/("d"k) = "d"/("d"k) int_0^(pi/2) dt/sqrt(1-k^2 sin^2 t)`
    `= int_0^(pi/2) k sin^2 t (1-k^2 sin^2 t)^(-3/2) dt`
    `= 1/k int_0^(pi/2) [(1-k^2 sin^2 t)^(-3/2) - (1-k^2 sin^2 t)^(-1/2)] dt`
    `= (E//{:k':}^2 - K)/k`.
  </span>
</p>

<p class="proposition">
  <b>Legendre 关系式</b>
  [菲赫金哥尔茨 第二册 511 目 12)]
  <span class="formula">
    `E K' + E' K - K K' = pi/2`.
  </span>
  其中 `K = K(k)`, `K' = K(k')`, `E = E(k)`, `E' = E(k')`.
</p>

<p class="proof">
  由
  <span class="formula">
    `("d"k')/("d"k) = "d"/("d"k) sqrt(1-k^2) = -k/sqrt(1-k^2) = -k/k'`,
  </span>
  将 `("d"E)/("d"k) = (E-K)/k`, `("d"K)/("d"k) = (E//{:k':}^2-K)/k`,
  `("d"E')/("d"k) = (-k/k')(E'-K')/k'`, `("d"K')/("d"k) = (-k/k')(E'//k^2 - K')/k'`
  代入整理得
  <span class="formula">
    `"d"/("d"k) (E K' + E' K - K K') = 0`.
  </span>
  即 `E K' + E' K - K K'` 等于一个常数 `c`. 现在令 `k to 0`, `k' to 1`
  求极限, 来求出 `c`. 事实上
  <span class="formula">
    `lim_(k to 0) E' K = 1 * pi/2 = pi/2`,
  </span>
  又
  <span class="formula">
    `|K'(E - K)|`
    `= |int_0^(pi/2) dt/sqrt(1-{:k':}^2 sin^2 t)| |int_0^(pi/2) (1/sqrt(1-k^2 sin^2 t) - sqrt(1-k^2 sin^2 t)) dt|`
    `le pi/(2 k) * |int_0^(pi/2) (k^2 sin^2 t dt)/sqrt(1-k^2 sin^2 t)|`
    `le pi/(2 k) * pi/2 * k^2/(k') to 0`,
  </span>
  所以极限为 `pi/2`.
</p>

<ol class="proposition">
  `|z| lt 1` 时, 椭圆积分展开为幂级数:
  <li>`K(z) = pi/2 sum_(n ge 0) [((2n-1)!!)/((2n)!!) z^n]^2`;</li>
  <li>`E(z) = -pi/2 sum_(n ge 0) 1/(2n-1) [((2n-1)!!)/((2n)!!) z^n]^2`.</li>
</ol>

<p class="proof">
  只证第一式. 利用幂级数
  <span class="formula">
    `1/sqrt(1-z) = sum_(n ge 0) ((2n-1)!!)/((2n)!!) z^n`,
    `quad |z| lt 1`
  </span>
  和积分
  <span class="formula">
    `int_0^(pi/2) sin^(2n) theta "d"theta = pi/2 ((2n-1)!!)/((2n)!!)`
  </span>
  对 `K(z)` 逐项积分即可.
</p>

<h3>第三类椭圆积分</h3>

<script src="../../js/note.js?type=math"></script>
<script src="../../js/plot.js"></script>
<script>
function agm (a, b) {
  while (Math.abs(a-b) > 1e-6) {
    var tmp = Math.sqrt(a*b)
    b = (a+b)/2
    a = tmp
  }
  return a
}
function ellipticK(k) {
  return Math.PI/(2 * agm(1, Math.sqrt(1-k*k)))
}
new Plot('plt-elliptic', {keepRatio: false, xmin: -1, xmax: 1, ymin: -0.5, ymax: 4}).axis().plot(ellipticK, {min: -1, max: 1})
</script>
</body>
</html>
